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Abstract. Optical flow cannot be computed locally, since only one independent measurement is available 
from the image sequence at a point, while the flow velocity has two components. A second constraint 
is needed. A method for finding the optical flow pattern is presented which assumes that the apparent 
velocity of the brightness pattern varies smoothly almost everywhere in die image. An iterative implemen- 
tation is shown which successfully computes the optical flow for a number of synthetic image sequences. 
The algorithm is robust in that it can handle image sequences that arc quantized rather coarsely in space 
and time. It is also insensitive to quantization of brightness levels and additive noise. Examples are 
included where the assumption of smoothness is violated at singular points or along lines in die image. 
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1. Introduction 

Optical flow is the distribution of apparent velocities of movement of brightness patterns in an image. 
Optical flow can arise from relative motion of objects and the viewer (Gibson 1950, 1966). Consequently, 
optical flow can give important information about the spatial arrangement of the objects viewed and the 
rate of change of this arrangement (Gibson 1977). Discontinuities in the optical flow can help in segment- 
ing images into regions that correspond to different objects (Nakayama & Loomis 1974). Attempts have 
been made to perform such segmentation using differences between succesive image frames (Jain et al. 
1977, Jain et al 1979, Jain & Nagcl 1979, Limb & Murphy 1975, Nagcl 1977). Some recent papers have 
considered the problem of recovering the motions of objects relative to the viewer from the optical flow 
(Hadani et al 1980, Kocndcrink & van Doom 1975 & 1976, Longuet-Higgins & Prazdny 1979, Prazdny 
1979 & 1980). In some cases information about the shape of an object may also be recovered (Kocndcrink 
& van Doom 1975 & 1976, Clocksin 1978). 

These papers begin by assuming that the optical flow has already been determined. Although some 
reference has been made to schemes for computing the fhw from successive views of a scene (Fcnncma 
& Thompson 1979, Hadani et al 1980), the specifics of a scheme for determining the flow from the image 
have not been described. Related work has been done in an attempt to formulate a model for the short 
range motion detection processes in human vision (Batali & Ullman 1980, Marr & Ullman 1979). 

The optical flow cannot be computed at a point in the image independently of neighboring points 
without introducing additional constraints, because the velocity field at each image point has two com- 
ponents while the change in image brightness at a point in the image plane due to motion yields only one 
constraint. Consider for example a patch of a pattern where brightness varies as a function of one image 
coordinate but. not the other. Movement of the pattern in one direction alters the brightness at a particular 
point, but motion in the other direction yields no change. Thus components of movement in the latter 
direction cannot be computed locally. Additional constraints must be introduced to fully determine the 
flow. 



2. Relationship to Object Motion 

The relationship between the optical flow in the image plane and the velocities of objects in the three 
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dimensional world is not necessarily obvious. We perceive motion when a changing picture is projected 
onto a stationary screen, for example. Conversely, a moving object may give rise to a constant brightness 
pattern. Consider for example, a uniform sphere which exhibits shading because its surface elements are 
oriented in many different directions. Yet, when it is rotated, the optical flow is zero at all points in the 
image, since the shading does not move with the surface. Also, specular reflections move with a velocity 
characteristic of the virtual image, not the surface in which light is reflected. 



3. The Restricted Problem Domain 

For convenience, we tackle a particularly simple world where the apparent velocity of brightness 
patterns can be directly identified with the movement of surfaces in the scene. To avoid variations in 
brightness due to shading effects we initially assume that the surface is flat. We further assume that 
die incident illumination is uniform across the surface. The brightness at a point in the image is then 
proportional to the reflectance of the surface at the corresponding point on the object. Also, we assume 
initially that reflectance varies smoothly and has no spatial discontinuities. This latter condition assures us 
that the image brightness is diffcrentiable. 

In this simple situation, the motion of the brightness patterns in the image is directly determined by 
the motions of corresponding points on the surface of the object. Computing the velocities of points on 
the object is a matter of simple geometry once the optical flow is known. 



4. Constraints 

We will derive an equation that relates the change in image brightness at a point to the motion of the 
brightness pattern. Let the image brightness at the point (x, y) in the image plane at time t be denoted by 
E(x, y, t). Now consider what happens when the pattern moves. The brightness of a particular point in 
the pattern is constant, so that 

dt 
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Using the chain rule for differentiation we see that, 

dEdx .dEdy dE 

dx dt + dy dt + dt ~~ ' 

(See Appendix A for a more detailed derivation.) If we let 

dx , dy 

u= Tt and „=-, 

then it is easy to see that we have a single linear equation in the two unknowns u and v, 

E x u + E y v+E t = Q, 

where we have also introduced the additional abbreviations E x , E y , and E t for the partial derivatives 
of image brightness with respect to x, y and t, respectively. The constraint on the local flow velocity 
expressed by this equation is illustrated in Figure 1. Writing the equation in still another way, 

{E x ,E y )-{u,v) = -E t . 

Thus the component of the optical flow in the direction of the brightness gradient {E x ,E y ) equals 



We cannot, however determine the component of movement in the direction of the iso-brightness con- 
tours, at right angles to the brightness gradient. As a consequence, the flow velocity (u, v) cannot be 
computed locally without introducing additional constraints. 



5. The Smoothness Constraint 

If every point of the brightness pattern can move independently, there is little hope of recovering the 
velocities. More commonly we view opaque objects of finite size undergoing rigid motion or deformation. 
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Figure 1. The basic rate of change of image brightness equation constrains the 
optical flow velocity. The velocity (u, v) has to lie along a certain line perpendicular 
to the brightness gradient vector (E x ,E y ) in velocity space. 
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In this case neighboring points on the objects have similar velocities and the velocity field of the bright- 
ness patterns in the image varies smoothly almost everywhere. Discontinuities in flow can be expected 
where one object occludes another. 

One way to express the additional constraint is to limit the difference between the flow velocity at 
a point and the average velocity over a small neighborhood containing the point. Equivalently we can 
minimize the sum of the squares of the Laplacians of the x- and y-components of the flow. The Laplacians 
of u and v are defined as 



„ 2 d 2 u , d 2 u J _ 2 d 2 v . &v 

Vu = £*+5y* and Vw = £*+V- 



In simple situations, both Laplacians are zero. If the viewer translates parallel to a flat object, rotates about 
a line perpendicular to the surface or travels orthogonally to the surface (assuming perspective projection), 
then the second partial derivatives of both u and v vanish. Note that our approach is in contrast with 
that of Fennema & Thompson (1979), who propose an algorithm that indirectly incorporates additional 
assumptions such as surface smoothness or object rigidity. 



6. Quantization and Noise 

Images may be sampled at intervals on a fixed grid of points. While tcssclations other than die 
obvious one have certain advantages (Mcrsercau 1979, Gray 1971), for convenience we will assume that 
the image is sampled on a square grid at regular intervals. Let the measured brightness be Eij t k at the 
intersection of the i-th row and jf-th column in the k-th image frame. Ideally, each measurement should 
be an average over die area of a picture cell and over the length of the time internal. In the experiments 
cited here we have taken samples at discrete points in space and time, instead. 

In addition to being quantized in space and time, the measurements will in practice be quantized in 
brightness as well. Further, noise will be apparent in measurements obtained in any real system. 
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7. Estimating the Partial Derivatives 

Wc must estimate the derivatives of brightness from the discrete set of image brightness measure- 
ments available. It is important that the estimates ofii^, E y , andE t be consistent. That is, they should 
all refer to the same point in the image at the same time. While there arc many formulas for approximate 
differentiation (Conte & de Boor 1965, Hamming 1962) wc will use a set which gives us an estimate of£ x , 
E y , E t at a point in the center of a cube formed by eight measurements. The relationship in space and 
time between these measurements is shown in Figure 2. Each of the estimates is the average of four first 
differences taken over adjacent measurements in the cube. 

E x ^-{Eij+i^ — Eij t k + J57*-|-i,iH-i,iB — Ei-\-i,j,k 

Ey ^tW'+i,^ ~ ^*,i»fc ~r" •£t+i,jH-i,fc ~~ Ei,j-\-i,k 

-f^'-f-ij^+i — E it j t k-\-Y H--E , *-fi,i+i,fc-+-i — ■Eij-f-i.fc+i} 

Et ^~{E i} j t k+\ —Ei,j,k -j-£t+i,i,JH-i — -£^-f-i t j,As 

"F'-^ij + l,*:+l — Ei,j-+-L,k H - -Ei-f-1 ,j-f-l,fc-h L — E{ -|-l,i+l,fc}- 

Here die unit of length is the grid spacing interval in each image frame and the unit of time is the image 
frame sampling period. 



8. Estimating the Laplacian of the Flow Velocities 

We also must approximate the Laplacians of u and v. One convenient approximation takes the 
following form 

V 2 u ra K,(u itj>k — Uij fk ) and V 2 v ^ K[v itjik — v it j ik ), 

where the local averages u and are defined as follows 

UiJ,k =g{Wt — lj,fc + Ui,j+\,k + «t-fl,i,fc + u i,j— l,k} 

+ yr{«»_l,j_ l,fc + v i— l,j+l,k + W»-|-lj+l,A: + v i-+-l,j— l,fc}- 
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Figure 2. The three partial derivatives of image brightness at the center of the cube 
are each estimated from the average of first differences along four parallel edges 
of the cube. Here the column index j corresponds to the x direction in the image, 
the row index i to the r/ direction, while k lies in the time direction. 
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The proportionality factor /c equals 3 if the average is computed as shown and assuming again that the 
unit of length equals the grid spacing interval. Figure 3 illustrates the assignment of weights to neighbor- 
ing points. 

9. Minimization 

The problem then is to minimize the sum of the errors in the equation for the rate of change of image 
brightness, 

S b = E x u+E y v+E t , 

and the estimate of the departure from smoothness in the velocity flow, 

& c = (u-u) 2 + (v-v) 2 . 

What should be the relative weight of these two factors? In practice the image brightness measurements 
will be corrupted by quantization error and noise so that we cannot expect &6 to be identically zero. This 
quantity will tend to have an error magnitude that is proportional to the noise in the measurement. This 
fact guides us in choosing a suitable weighting factor, denoted by a 2 , as will be seen later. 
Let the total error to be minimized be 

g2 = , Q 2g2_ [ _g2 

The minimization is to be accomplished by finding suitable values for the optical flow velocity (u, v). We 
differentiate 8 2 to obtain 

,9£ 2 

V- = - 2a 2 (u - u) + 2{E x u + E y v + E t )E x 

ou 

^- = - 2a 2 (t) - v) + 2{E x u -f E y v + E t )E y . 
ov 

Setting these two derivatives equal to zero leads to two equations in u and v, 

{a 2 + E 2 x )u + E x E u v = (o 2 u - E x E t ) 
EJSyU + (a 2 + E 2 y )v = (a 2 i> - EyE t ). 
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Figure 3. The Laplacian is estimated by subtracting the value at a point from a 
weighted average of the values at neighboring points. Shown here are suitable 
weights by which values can be multiplied. 
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The determinant of the coefficient matrix equals a 2 {a 2 -\-E 2 x + E y ). Solving for u and v we find that 

(a 2 + E\ + E\ )u = + (a 2 + £ 2 )u - E^ - £*£, 
(a 2 + £ 2 + E 2 y )v = - E x E y u + (a 2 + ^)t) - EyE t . 



10. Difference of Flow at a Point from Local Average 

These equations can be written in the alternate form 

(a 2 + E\ + JE*)(u - u) = - £*[£,& + E y v + £ t ] 
(a 2 + El + E 2 )(<; - 0) = - £^ti + £ y t> -f- £]. 

This shows that the value of the flow velocity (u, v) which minimizes the error 6 2 lies in the direction 
towards die constraint line along a line that intersects die constraint line at right angles. This relationship 
is illustrated geometrically in Figure 4. The distance from die local average is proportional to die error 
in die basic formula for rate of change of brightness when it, v arc substituted for u and v. Finally we 
can see that a 2 plays a significant role only for areas where the brightness gradient is small, preventing 
haphazard adjustments to die estimated flow velocity occasioned by noise in the estimated derivatives. 
This parameter should be roughly equal to die expected noise in the estimate of E 2 -f E 2 . 

11. Constrained Minimization 

When we allow a 2 to tend to zero we obtain the solution to a constrained minimization problem. 
Applying the method of Lagrange multipliers (Russell 1976, Yourgau & Mandclstam 1968) to the 
problem of minimizing g 2 . while maintaining 6^ = leads to 

(£* "f E y ){u -u) = - E x [E x u + Eyv + E t \ 
{E\ + E 2 y )(v -v) = - E y [E x u + E y v + E t \. 

We will not use dicse equations since we do expect errors in the estimation of the gradient (E X) E y ). 
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Figure 4. The value of the flow velocity which minimizes the error lies on a line 
drawn from the local average of the flow velocity perpendicular to the constraint 
line. 
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12. Iterative Solution 

We now have a pair of equations for each point in the image. It would be very costly to solve these 
equations simultaneously by one of the standard methods, such as Gauss- Jordan elimination (Hamming 
1962, Hildebrand 1952). The corresponding matrix is sparse and very large since the number of rows and 
columns equals twice the number of picture cells in the image. Iterative methods, such as the Gauss-Seidel 
method (Hamming 1962, Hildebrand 1956), suggest themselves. We can compute a new set of velocity 
estimates (tt""* -1 , i> n_t_1 ) from the estimated derivatives and the average of the previous velocity estimates 
[u n , v n ) by 

u »+l = fc» _ Ex[Exil n + £ y - „ + Ei]/{q 2 + E 2 + E 2 y) 

„»+i =0" - E y \E x u n + Eyv n + E t ]/(a 2 + E 2 X + E*). 

It is interesting that the new estimates at a particular point do not depend directly on the previous 
estimates at the same point. 



13. Filling In Uniform Regions 

In parts of the image where the brightness gradient is zero, the velocity estimates will simply be 
averages of the neighboring velocity estimates. There is no local information to constrain the apparent 
velocity of motion of the brightness pattern in these areas. Eventually the values around such a region will 
propagate inwards. If the velocities on the border of the region arc all the same, dien points in the region 
will be assigned that value too after a sufficient number of iterations. Velocity information is thus filled in 
from the boundary of a region of constant brightness. 



14. Number of Iterations 

If the values on the border are not all the same, the values filled in will correspond to the solution 
of die Laplace equation for the given boundary condition (Ames 1977, Milne 1953, Richtmycr& Morton 
1957). The progress of this filling in phenomena is similar to the propagation effects in the solution of (lie 
heat equation, where the time rate of change of temperature is proportional to the I.aplacian. This gives 
us a means of understanding the iterative method in physical terms and of estimating the number of steps 
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required. The number of iterations should be larger than the cross-section of the biggest region that must 
be filled in. If the size of such regions is not known in advance one may use the cross-section of the whole 
image as a conservative estimate. 



15. Tightness of Constraint 

When brightness in a region is a linear function of the image coordinates we can only obtain the 
component of optical flow in the direction of the gradient. The component at right angles is filled in from 
the boundary of the region as described before. In general the solution is most accurately determined 
in regions where the brightness gradient is not too small and varies in direction from point to point. 
Information which constrains both components of the optical flow velocity is then available in a relatively 
small neighborhood of each point. Too violent fluctuations in brightness on the other hand are not 
desirable since the estimates of derivatives will be corrupted due to undcrsampling and aliasing. 



16. Choice of Iterative Scheme 



As a practical matter one has a choice of how the iterations are to be interlaced with the time steps. 
On the one hand, one could iterate until the solution has stabilized before advancing to the next image 
frame. On die other hand, given a good initial guess one may need only one iteration per time-step. A 
good initial guess for the optical flow velocities is usually available from the previous time-step. 

The advantages of the latter approach include an ability to deal with more images per unit time and 
better estimates of optical flow velocities in certain regions. Areas in which die brightness gradient is small 
lead to uncertain, noisy estimates obtained partly by filling in from the surroundings. These estimates arc 
improved by considering further images. The noise in measurements of the images will be independent 
and tend to cancel out. Perhaps more importantly, different parts of the pattern will drift by a given point 
in the image 

A practical implementation would most likely employ one iteration per time step for these reasons. 
We illustrate both approaches in the experiments. 
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17. Experiments 

The iterative scheme has been implemented and applied to image sequences corresponding to a 
number of simple flow patterns. The results shown here are on a relatively small image of 32 by 32 
picture cells. The brightness measurements were intentionally corrupted by approximately 1% noise and 
then quantized into 256 levels to simulate a real imaging situation. The underlying surface reflectance 
pattern was a linear combination of spatially orthogonal sinusoids. Their wavelength was chosen to give 
reasonably strong brightness gradients without leading to undcrsampling problems. Discontinuities were 
avoided to ensure that the required derivatives exist everywhere. 

Shown in Figure 5 for example are four frames of a sequence of images depicting a sphere rotating 
about an axis inclined towards the viewer. A smoothly varying reflectance pattern is painted on the surface 
of die sphere. The sphere is illuminated uniformly from all directions so that there is no shading. 



18. Results 

The first flow to be investigated was a simple linear translation of the entire brightness pattern. The 
resulting computed flow is shown as a needle diagram in Figure 6 for 1, 4, 16, and 64 iterations. The 
estimated flow velocities are depicted as short lines, showing the apparent displacement during one time 
step. In this example a single time step was taken so that the computations are based on just two images. 
Initially the estimates of flow velocity are zero. Consequcndy the first iteration shows vectors in die 
direction of the brightness gradient. Later, the estimates approach the correct values in all parts of die 
image. Few changes occur after 32 iterations when the velocity vectors have errors of about 10%. The 
estimates tend to be too small, rather than too large, perhaps because of a tendency to underestimate the 
derivatives. The worst errors occur, as one might expect, where die brightness gradient is small. 

In the second experiment one iteration was used per time step on the same linear translation image 
sequence. The resulting computed How is shown in Figure 7 for 1, 4, 16, and 64 time steps. The estimates 
approach die correct values more rapidly and do not have a tendency to be too small, as in the previous 
experiment. Few changes occur after 16 iterations when die velocity vectors have errors of about 7%. The 
worst errors occur, as one might expect, where the noise in recent measurements of brightness was worst. 
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Figure 5. Four frames out of a sequence of images of a sphere rotating about an 
axis inclined towards the viewer. The sphere is covered with a pattern which 
varies smoothly from place to place. The sphere is portrayed against a fixed, lightly 
textured background. Image sequences like these are processed by the optical flow 
algorithm. 
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Figure 6. Flow pattern computed for simple translation of a brightness pattern. The 
estimates after 1, 4, 16, and 64 iterations are shown. The velocity is 0.5 picture 
cells in the x direction and 1.0 picture cells in the y direction per time interval. 
Two images are used as input, depicting the situation at two times separated by 
one time interval. 
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Figure 7. Flow pattern computed for simple translation of a brightness pattern. The 
estimates after 1, 4, 16, and 64 time steps are shown. Here one iteration is used 
per time step. Convergence is more rapid and the velocities are estimated more 
accurately. 
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While individual estimates of velocity may not be very accurate, the average over the whole image was 
within 1% of the correct value. 

Next, the method was applied to simple rotation and simple contraction of the brightness pattern. 
The results after 32 time steps are shown in Figure 8. Note that the magnitude of the velocity is propor- 
tional to the distance from the origin of the flow in both these cases. 

In all of the examples so far the Laplacian of both flow velocity components is zero everywhere. We 
also studied more difficult cases where this was not the case. In particular, if we let the magnitude of the 
velocity vary as the inverse of the distance from the origin we generate flow around a line vortex and two 
dimensional flow into a sink. The computed flow patterns are shown in Figure 9. Here the computation 
involved many iterations based on a single time step. The worst errors occur near die singularity at the 
origin of the flow pattern, where velocities are found which arc much larger than one picture cell per time 
step. 

Finally we considered rigid body motions. Shown in Figure 10 are the flows computed for a cylinder 
rotating about its axis and for a rotating sphere. In both cases the Laplacian of the flow is not zero and in 
fact the Laplacian of one of the velocity components becomes infinite on the occluding bound. Since the 
velocities themselves remain finite, reasonable solutions are still obtained. The correct flow patterns are 
shown in Figure 11. Comparing the computed and exact values shows that the worst errors occur on die 
occluding boundary. These boundaries constitute a one dimensional subset of the plane and so one can 
expect tiiat the relative number of points at which die estimated flow is seriously in error will decrease as 
the resolution of the image is made finer. 

In Appendix B it is shown that there is direct relationship between the Laplacian of die flow velocity 
components and die Laplacian of the surface height. This can be used to see how our smoothness con- 
straint will fair for different objects. For example, a rotating polyhedron will give rise to a flow which has 
zero Laplacian except on the image lines which arc die projections of die edges of die body. 



19. Summary 

A method was developed for computing optical flow from a sequence of images. It is based on the 
observation that the flow velocity has two components and that die basic equation for the rate of change 
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Figure 8. Flow patterns computed for simple rotation and simple contraction of a 
brightness pattern. In the first case, the pattern is rotated about 2.8 degrees per 
time step, while it is contracted about 5% per time step in the second case. The 
estimates after 32 time steps are shown. 
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Figure 9. Flow patterns computed for flow around a line vortex and two dimensional 
flow into a sink. In each case the estimates after 32 iterations are shown. 
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Figure 10. Flow patterns computed for a cylinder rotating about its axis and for a 
rotating sphere. The axis of the cylinder is inclined 30 degrees towards the viewer 
and that of the sphere 45 degrees. Both are rotating at about 5 degrees per time 
step. The estimates shown are obtained after 32 time steps. 
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Figure 11. Exact flow patterns for the cylinder and the sphere. 
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of image brightness provides only one constraint. Smoothness of the flow was introduced as a second 
constraint. An iterative method for solving the resulting equation was then developed. A simple im- 
plementation provided visual confirmation of convergence of the solution in the form of needle diagrams. 
Examples of several different types of optical flow patterns were studied. These included cases where the 
Laplacian of the flow was zero and cases where it became infinite at singular points or along bounding 
curves. 

20. Acknowledgement 

One of the authors (Horn) would like to thank Professor H.-H. Nagel for his hospitality. The basic 
equations were conceived during a visit to the University of Hamburg, stimulated by Professor Nagel's 
long-standing interest in motion vision. The other author (Schunck) would like to thank W. E. L. Grimson 
and E. Hildreth for many interesting discussions and much knowledgable criticism. 



Horn & Schunck A. I. Memo 572 



Page 24 Determining Optical Flow 

Appendix A. Rate of Change of Image Brightness 

Consider a patch of the brightness pattern that is displaced a distance Sx in the ^-direction and Sy in 
the y-direction in time St. The brightness of the patch is assumed to remain constant so that 

E{x, y, t) = E{x + 6x,y + Sy, t -f St). 

Expanding the right hand side about the point (x, y, t) we get, 

E(x, y, t) = E(x, y, t) + 6x d £ + 5yg + St^ + e. 

Where e contains second and higher order terms in Sx, 6y, and St. After subtracting E{x, y, t) from both 
sides and dividing through by St we have 

SI +£f +!+««>-* 

where 0(St) is a term of order St, and we assume that Sx and Sy vary as St. In die limit as St -> this 
becomes 

dx dt + dy dt^ ' dt 
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Appendix B. Smoothness of Flow for Rigid Body Motions 



Let a rigid body rotate about an axis (u) x , uj y , uo z ), where the magnitude of the vector equals the 
angular velocity of the motion. If this axis passes through the origin, then the velocity of a point (x, y, z) 
equals the cross product of K, w y , u> z ), and (x, y, z). There is a direct relationship between the image 
coordinates and the x and y coordinates here if we assume that the image is generated by orthographic 
projection. The x and y components of the velocity can be written, 

u =u)yZ — uj z y 

V =Ij) z X — U) X Z. 

Consequently, 

V 2 u = + w y V 2 z 
V 2 v = — u) x V 2 z. 

This illustrates that the smoothness of the optical flow is related directly to the smoothness of the rotating 
body and that the Laplacian of the flow velocity will become infinite on the occluding bound, since the 
partial derivatives of* with respect to x and y become infinite there. 
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